Some remarks on the inverse Smoluchowski problem for cluster-cluster aggregation 
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It is proposed to revisit the inverse problem associated witli Smoluchowski's coagulation equation. 
The objective is to reconstruct the functional form of the collision kernel from observations of the 
time evolution of the cluster size distribution. A regularised least squares method originally proposed 
by Wright and Ramkrishna (1992) based on the assumption of self-similarity is implemented and 
tested on numerical data generated for a range of different collision kernels. This method expands 
the collision kernel as a sum of orthogonal polynomials and works best when the kernel can be 
expressed exactly in terms of these polynomials. It is shown that plotting an "L-curve" can provide 
an a-priori understanding of the optimal value of the regularisation parameter and the reliability of 
the inversion procedure. For kernels which are not exactly expressible in terms of the orthogonal 
polynomials it is found empirically that the performance of the method can be enhanced by choosing 
a more complex regularisation function. 



INTRODUCTION 

The effects of clouds and precipitation is one of the 
largest sources of uncertainty in our current attempts 
to simulate the Earth's climate [ij. The reason for this 
is that most phenomena associated with clouds happen 
on scales below those which are explicitly resolved by 
the current climate models. Their feedback on resolved 
scales must therefore be parameterised. Such parameter- 
isations require a strong understanding of the underlying 
physical processes taking place. One such process which 
has attracted considerable attention in recent years is the 
time evolution of the droplet size distribution in clouds 
and its connection with precipitation. Water droplets in 
clouds are seeded by cloud condensation nuclei! such as 
aerosol particles. They initially grow by condensation 
and subsequently by coagulation as droplets collide with 
each other in the cloud to produce larger droplets. The 
detailed micro-physics of the coagulation process is diffi- 
cult to understand theoretically since turbulence within 
the cloud is believed to play a key role in determining 
the collision rate of droplets A complete theoreti- 
cal description is hampered by the difficulty in describ- 
ing the statistical interplay between particles and tur- 
bulence analytically. On the other hand the quality of 
the data available for the study of this problem is im- 
proving rapidly due to recent advances in observational 
techniques Q and direct numerical simulation of droplets 
in turbulent flows 0, 0] . 

In the context of the droplet coagulation problem, 
much research has focused on using our improved ob- 
servational understanding of the behavious of droplets 
in turbulent flows to calculate more accurately the func- 
tional form of the droplet coagulation rate, X (mi, 7712), 
as a function of droplet masses, or equivalently (if 
droplets are assumed spherical), droplet size. In this 
note, we argue that the improved availability of data 
suggests that we should also revisit the corresponding 



inverse problem which can be stated as follows: given 
observations or measurements of the time evolution of 
the droplet size distribution, how much information can 
be extracted about the functional form of the coagulation 
rate? This problem was studied in the past by Wright 
and Ramkrishna 6] in the context of chemical mixing and 
has been recently revisited in the context of droplets in 
turbulence by Onishi et al. Q . Although there are many 
difficulties associated with such inverse problems, as dis- 
cussed below, there are potential rewards. The approach 
is data driven and could provide a guide to modeling in 
situations where the underlying microphysics remains in- 
completely understood. The inverse approach could also 
begin to quantify the extent to which available data and 
measurements can distinguish between different models. 

First a word about the usual forward problem. If the 
collision kernel, K{mi,m2), is known, the evolution of 
the cluster size distribution, Nm{t) is described by the 
Smoluchowski coagulation equation Q: 



dtN„i {t) = / dmi K{mi,m - mi) N„n (t) N„ 



2N„,{t) / dmiK{m,mi)N,n^{t). 



(1) 



It is applicable when spatial correlations between parti- 
cles are sufficiently weak that collisions between parti- 
cles can be considered statistically independent. A huge 
amount is known theoreticallyabout the solutions of the 
Smoluchowski equation. See [9j for a modern review. 



HOMOGENEOUS COLLISION KERNELS AND 
SELF-SIMILARITY 



In some applications (see [lO| for some discussions), 
the collision kernel is a homogeneous function of its ar- 
guments whose degree we shall denote by 7: 

K{ami, am2) = a? K{mi, 7712). 
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FIG. 1. Snapshots of the self-similar evolution of the clus- 
ter size distribution, Nm{t), obtained from a numerical sim- 
ulation of Eq. lU with the collision kernel K(mi,ra2) = 

(mi, 7712)1. The inset shows the scaling function obtained 
when the data is collapsed according to Eq. ^ with the typ- 
ical scale, s(t), obtained as in Eq. Q with n = 2. 



For such kernels, the evolution of the cluster size distri- 
bution is often self-similar. That is to say we can write: 



N^{t) ^ s{t)-^ f [z) 



m 



(2) 



where s{t) is the typical cluster size. This can be defined 
intrinsically as a ratio of moments of the size distribution: 



sit) 



Mn{t) 



Mp{t) 



mPNm(t)dm. (3) 



The scaling function, f{z), determines the shape of the 
cluster size distribution. An example of this self-similar 
evolution obtained from numerical simulation of Eq. ([!]) 
with K{mi,m2) = (mim2)^/'* is shown in Fig. [T] The 
inset shows the scaling function, f{z), obtained by col- 
lapsing the data according to Eq. ^ with the typical 
scale, s{t), obtained as in Eq. ^ with n = 2. All nu- 
merical simulations of Eq. ([1]) reported in this paper were 
done using the pairwise binning method described in [Tlj . 

A homogeneous collision kernel and self-similar time 
evolution are not plausible assumptions for the droplet 
coagulation problem in clouds which has motivated this 
study If one restricts attention to the gravitational 
coagulation-dominated regime, one could perhaps argue 
that the kernel is approximately homogeneous. It turns 
out, however, that even with the resulting differential 
sedimentation kernel which is homogeneous of degree 
7 = 4/3, the Smoluchowski equation does not produce a 
self-similar evolution but rather undergoes instantaneous 
gelation [l3, 14 1. Our focus on self-similar problems sim- 
ply stems from the fact that they provide the simplest 
class of coagulation problems on which we can begin to 
study the inverse problem described above. We further 



restrict ourselves to kernels having 7 < 1 in order to 
avoid dealing with the gelation transition. 



THE INVERSE SMOLUCHOWSKI PROBLEM 

Let us introduce the cumulative cluster mass distribu- 
tion: 



mi Njm(t)dmi. 



(4) 



The original cluster size distribution, Nm{t) is recovered 
from J-m{t) by differentiation: 



m om 



(5) 



In terms of Eq. ([T]) can be written in the compact 

form: 



1712 



K{mi, 771,2) 



(6) 

If we assume scahng as in Eq. ^ then Fm{t) = F{z) and 
the scaling function of the cumulative mass distribution, 
F(z), satisfies the following scaling equation: 



dF{z2) 

Z2 



K{zi,Z2). (7) 



This is a complicated nonlinear integro-differential equa- 
tion if we are required to determine F{z) given K{zi^ Z2)) 
(the forward problem). It is, however, a linear equation 
if we are required to determine K{zi, Z2) given F{z) (the 
inverse problem). This inverse problem is, however, ill- 
posed. This is most easily seen by considering what hap- 
pens if we discretise Eq. ^ on N z-points. We obtain a 
linear system, 

g = 5k, 

consisting of equations for the A^^ values of Ar(2i, Z2) 
on the discretisation points. Actually we can reduce 
the number of unknowns almost by a factor of 2 since 
K{zi,Z2) ~ A'(z2,zi) but the conclusion remains un- 
changed. This system is enormously under-determined, 
which implies that one can find many solutions but they 
are typically over-fitting the measurements of F{z). 

One way of dealing with under-determined systems 
is via Tikhonov Regularisation (Ridge regression). The 
idea is to solve a minimization problem with a regulari- 
sation term: 



argmm {HS'k-gll 



MM'} 



(8) 



Noise-dominated solutions have to compete against the 
regularization term A||k|p in the minimisation. This ap- 
proach was pioneered by Wright and Ramkrishna Q. It 
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has the advantage that, in principle, one does not need to 
know the functional form of the collision kernel a-priori. 
The trick is to choose the "best" value of the regulariza- 
tion parameter, A. Before discussing the selection of the 
value of A, we discuss briefly how to set up the minimiza- 
tion problem. 

Following 0, we assume that K {mi, -012) can be ex- 
panded in terms of Laguerre polynomials, Li(x), up to 
order p in the mi and directions: 

if(mi,m2) = ^a„L„(mi,m2) (9) 

n=l 

where 

Ln{x,y) = Li{x)Lj{y) 

where n = {i — 1) p + j is a compound index which rolls 
the indices, i and j, in the mi and m2 directions into 
one. With Eq. ([9]), the minimisation problem Eq. dH) is 
transformed into a different minimisation which aims to 
determine the coefficients, a„ in Eq. ([9]). This has a num- 
ber of advantages. It ensures that the minimisation prob- 
lem automatically finds kernels which are fairly smooth 
functions of mi and m2. It also enormously reduces the 
size of the problem since the number of polynomials, p, 
required in each direction is typically rather small. Fur- 
thermore, this choice allows the entries of the matrix S in 
Eq. dS]) to be calculated semi-analytically as described in 
0. It has the disadvantage, however, of constraining the 
solution to those kernels which can be expressed in the 
form of Eq. (|9]) which effectively requires a-priori knowl- 
edge of the class of kernels which we are seeking, as was 
done in Q. 

As mentioned above, the tricky part of this procedure 
is to determine the most apprpriate value of A to use 
in Eq. (jS]). If A is too small, the result will over- fit the 
observations. If A is too large then kest will be pushed 
to zero by the weight of the second term in Eq. ^ and 
will retain little information about the Smoluchowski dy- 
namics encoded in the matrix S. A rational approach to 
determining A is provided by plotting an "L-curve" (for 
a clear review see This is a plot of the norm of the 

solution, ||k||, as a function of the norm of the residual, 
||Ak — b||, for a range of values of A. If there is a clear 
transition in the minimisation problem from a regime 
dominated by overfitting to a regime dominated by the 
rcgularisation then the L-curve will have the distinctive 
L-shape and the best values of A are those in the "elbow" 
of the curve. 

RESULTS 

The regularised least squares method described above 
was applied to several sets of data obtained by the nu- 
merical integration of Eq. ([T]) with monodisperse initial 



data with different model kernels. We used 40 discretisa- 
tion points in the interval z G [0.01, 2.0] and chose p = 3 
Laguerre polynomials in each direction. In each case, the 
data collapse and extraction of the scaling function was 
done using the method described in [l^ and then fitted 
to a function of the form Cz"~^e~^^ as suggested in 
0. These fitted scaling functions were then used as the 
"observations", F{z), in the discretisation of Eq. ([7]). 

Figs. [2] and |3] show some results of the regularised least 
squares for the constant kernel, K(mi,m2) = 1, and 
the sum kernel, K{mi,m2) = ^(mi + m2), respectively. 
These are good test cases to begin with since the solu- 
tion of Eq. ([T]) with monodisperse initial data is known 
explicitly for each of these kernels allowing the numeri- 
cal aspects of the calculation to be validated. The right 
panels of the figures show the L-curves obtained by per- 
forming the minimisation (jS]) over a range of values of 
A. In both cases, a clear elbow is visible in the result- 
ing curve indicating an appropriate range of values for 
A. The left panel shows, for ease of visualisation, the di- 
agonal, K{m,m), of the kernels obtained by regularised 
least squares as a function of m for different values of A. 
The A curves shown correspond to values of A which are 
too big, too small and "just right" (meaning a value in 
the elbow of the L-curve). It is clear that the method 
does a reasonable job of extracting the basic shape of the 
kernel given the scaling function, F(z). Surprisingly, we 
found that it was not necessary to explicitly enforce the 
positivity of K{mi^2) by constrainting the minimisa- 
tion as was done in [6[ . The method seemed to do equally 
well, and in some cases, better without these constraints. 
For the sum kernel, we did enforce the constraint that 
the kernel should vanish for particles of zero mass which 
seems physically reasonable - the results were less con- 
vincing without this constraint. 

Fig. 0] shows the corresponding results for the gener- 
alised sum kernel K{mi,m2) = ^{y/mi + ^m2 ) . The 
results are clearly less convincing. This is probably re- 
lated to the fact that this kernel cannot be expressed in 
terms of the Laguerre polynomials chosen to represent 
the kernel in Eq. ©. It is worth pointing out that the 
corresponding L-curve does not have a sharp elbow indi- 
cating that the method does not find a clear "best" value 
of A. This seems to illustrate a point in favour of such 
methods - the fact that the L-curve does not have a sharp 
elbow provides an a-priori indication that the results of 
the minimisation should be treated with caution. 

Empirical investigation suggest that the results for the 
generalised sum kernel can be improved by tinkering with 
the rcgularisation procedure. If, instead of, Eq. we 
perform the minimisation 

mm\\Sk-g\\l + Xw{k) (10) 
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FIG. 2. Results for the constant kernel, K{mi,m2 = 1. The right panel shows the L-curve obtained by performing the 
minimisation Eq. ([8]) over a range of values of A. The left panel shows the diagonal, K{z,z), of the reconstructed kernels 
compared to the theoretical curve for values of A in the upper left, lower right and the "elbow" of the L-curve. 
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FIG. 3. Results for the sum kernel, K{mi,ni2 = \{m\ + 7712). The right panel shows the L-curve obtained by performing 
the minimisation Eq. ([8]) over a range of values of A. The left panel shows the diagonal, K{z, z), of the reconstructed kernels 
compared to the theoretical curve for values of A in the upper left, lower right and the "elbow" of the L-curve. 




FIG. 4. Results for the generalised sum kernel, K{mi,m2 = ^(y^mT+ ^m,2). The right panel shows the L-curve obtained by 
performing the minimisation Eq. ((8| over a range of values of A. The left panel shows the diagonal, K{z, z), of the reconstructed 
kernels compared to the theoretical curve for values of A in the upper left, lower right and the "elbow" of the L-curve. 
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where 



;(k)=^log(|fc. 



1) = loe 



n( 



1). 



(11) 



wc found that the results were much better. This form 
for w was found by experimentation. At present, this 
statement is at the level of empirical observation and 
requires further investigations. 

DISCUSSION AND OUTLOOK 

Wc conclude, as several previous authors have done, 
that the inverse Smoluchowski problem is technically 
feasible and could potentially be developed into a use- 
ful tool for the study of droplet size distributions in 
clouds and other applications where the underlying mi- 
crophysics is still incompletely understood. The results 
presented here indicate that the L-curve provides a useful 
complementary tool to the methods developed by Wright 
and Ramkrishna @ to allow the regularisation parame- 
ter to be selected a-priori in situations where the collision 
kernel is not known from the outset. It is worth men- 
tioning that the results presented in Figs. ©-Q do not 
give a good indication to the casual reader of the degree 
of numerical sensitivity required in tackling these prob- 
lems. It became clear to us during these investigations 
that the ill-posedness of the inverse problem represented 
by Eq. ([7]) requires that great care be taken in interpret- 
ing the outputs of these methods. It is also clear that 
further research is required, even in the simple case of 
self-similar time evolution, in order to make the method 
more robust. Our results suggest that we should consider 
more general functional forms for the collision kernel than 
Eq. ^ in order to improve this robustness. This will 
probably not pose much difficulty since the analytic sim- 
plification obtained by the use of Laguerre polynomials 
in is probably less important nowadays owing to the 
increased computational power which can be brought to 
bear on the computation of matrix elements by quadra- 
ture when closed analytic forms are not available. 

In the long run, however, it is clear that it is necessary 
to free ourselves of the assumptions of homogeneity of 



the kernel and self-similarity of the cluster size distribu- 
tion. Some strong progress in this direction has already 
been made recently by Onishi et al. who have been 
able to infer the relative importance of the turbulent and 
gravitational coagulation as a function of Reynolds num- 
ber from direct numerical simulations of droplet-laden 
turbulence using inverse methods. Our approach differs 
slightly from this work in the sense that we would like, as 
far as possible, to learn the functional shape of the ker- 
nel from the data by allowing considerable freedom in the 
class of possible kernel functions. This approach would 
be more appropriate in situations when the underlying 
micro-physics is unknown or controversial. 
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